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Abstract 


The truncation theory as it pertains to the calculation of geo id undulations 
based on Stokes’ integral, but from limited gravity data, is reexamined. Speci- 
fically, the improved procedures of Molodenskii et al, (1962) are shovyrn through 
numerical investigations to yield substantially smaller errors than the conven- 
tional method that is often applied in practice. In this improved method, as well 
as in a simpler alternative (Meissl, 1971b) to the conventional approach, the 
Stokes' kernel is suitably modified in order to accelerate the rate of convergence 
of the error series. These modified methods, however, effect a reduction in the 
error only if a set of low -degree potential harmonic coefficients is utilized in 
the computation. Consider, for example, the situation in which gravity anomalies 
are given in a cap of radius 10° and the GEM 9 (20,20) potential field is used. 
Then, typically, the error in the computed undulation (aside from the spherical 
approximation and errors in the gravity anomaly data) according to the conven- 
tional truncation theory is 1.09 m; with Meissl’s modification it reduces to 
0.41 m, while Molodenskii's improved method gives 0.46 m. A further altera- 
tion of Molodenskii's mediod is developed and yields an RMS error of 0. 33 m. 
These values reflect the effect of the truncation, as well as the errors in the 
GEM 9 harmonic coefficients. The considerable improvement, suggested by 
these results, of the modified methods over the conventional procedure is verified 
with actual gravity anomaly data in two oceanic regions, where the GEOS -3 altim- 
eter geoid serves as the basis for comparison. The optimal method of truncation, 
investigated by Colombo (1977), Is extremely ill-conditioned. It is shown ihat 
with no corresponding regularization, this procedure is inapplicable. 
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1. Introductton 


The practical application of Stokes' solution to the geodetic boundary*value 
problem (Stokes' integral formula) has been studied extensively. The principal 
difficulty with this solution in practice is the formal requirement of continuous 
gravity data covering the entire geoid (approximated by a sphere). The lack of 
global coverage has led to the natural approximation whereby the integration is 
limited to a spherical cap ("truncation" of the integral). However» the data within 
a cap do not always constitute the total input, since much of the gravity information 
of the remaining exterior areas is often conveniently derived from a number of 
low-degree (say, up to degree 20) harmonic constituents which are determined, 
for example, by satellite methods. 

The theory of this combination of satellite and terrestrial data to obtain the 
disturbing potential (or geoid undulation) originates with M. S. Molodenskii (Molo- 
denskii, 1958; Molodenskii et al. , 1962). The derivation of Molodenskii's basic 
theory has also been carried out by Heiskanen and Moritz (1967), and their expo- 
sition is the most familiar. However, Molodenskii (1958) (see also Molodenskii 
et al. , 1962) extended and refined his theory leading to significant reductions in 
the truncation error (presuming a knowledge of a set of potential harmonic coeffi- 
cients). In the recent literature, several authors have become aware of potential 
reductions in the error (Meissl, 1971b; Dickson, 1979), but Molodenskii must be 
credited with the most rigorous and comprehensive treatment. The applications of 
this theory, including further possible improvements, have been studied by Hsu 
Houtze (Hsu Houtze and Zhu Zhuowen, 1979). 

It is the purpose of this report to bring to light Molodenskii’s theory which 
seems to have gone largely unnoticed by practicing geodesists and to test its appli- 
cability with the current data available. In this respect, the gravity data inside the 
cap are assumed to be error free and sufficiently dense so that numerical integration 
errors can be neglected; on the other hand, errors in the harmonic coefficients of 
the gravity potential are allowed and treated accordingly. 

Meissl (1971b) has suggested a very simple modification (which has also been 
used by Ostach (1970), see also Hsu Houtze and Zhu Zhuowen, 1979) to Stokes' func- 
tion, effecting a substantial reduction in the truncation error if a higher-degree 
reference field is given. Both Meissl's approach and Molodenskii's original ideas 
rely on identical principles, but Molodenskii specifically builds on the premise that 
the geoid undulation computation should utilize potential harmonic coefficients. 

Meissl does not mention Ihe use of harmonic coefficients and in applications his 
modification is quite useless in their absence. Nevertheless, since Meissl presents 
considerable mathematical insight into the truncation of Stokes' integral, his treat- 
ment will also be discussed in detail. The investigations would not be complete 
without mentioning the work of Wong and Gore (1969) whose method of introducing 
harmonic components is, however, somewhs;t arbitrary. Further tests of their 
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method are conducted in this report. Finals , the entirely different aj^roach to 
die truncation theory of Colombo <1977) merits investigation. Colombo's solution, 
although yielding promising results, is only approximate because of the numerical 
instability of the problem posed by him. A rigorous solution is presented here, 
which however does not bypass the instability. 

Most of the relevant derivations of Molodenskii, Melssl, and Wong and Gore 
ai« reproduced In the following sections, the only intention being to make this report, 
to some degree, self-contained. It is attempted to maintain a simple notation. Un- 
fortunately, since a variety of kernels and modifications of kernels will enter the 
discussion, the symbolism must deviate occasionally from the usual notation found 
in the literature. The following will be adhered to; 

Ag geoidal gravity anomaly 

Ag„ nth surface harmonic function of Ag 

estimate of Ag„ based on known liarmonic coefficients (which may 
be subject to error) 

6(Ag„) the error in Ag„ such that Ag„ = Afe„ + 6<Ag„) 
c„ degree variance of Ag„ 

6c n degree variance of 6(Agn) 

N geo id undulation given by Stokes' integral 

1^ undulation computed by integrating over a cap, and possibly using 
a number of harmonic coefficients 

6n the error in ^I, such that N = + 6 n 

3n the RMS (root mean square, global average) of 6N 

W the RMS error, specifically when a number of harmonic coefficients 

are known 

S(y) Stokes' function, y = cos <|!) 

S*(y) S(y) minus its harmonics up to degree m 

S^(y) the first h -f 1 harmonics of S when the latter is expanded in the 
interval [-l,yol, yo< 1 

K(y) the kernel used for the integration over the cap 

AK(y) the error kernel of the integral representing 6N 

Q„ a truncation coefficient (for any one of the kernels) 

The individual modifications to the "classic” truncation theory are represented by 
subscripts, the subscript 1 referring to the classic theory. A subscript 2 denotes 
Meissl's (1971b) modification, a 3 refers to the method of Wong and Gore (1969), 
and an M signifies Molodenskil's (1958) modification. Other notations are explained 
in the text. 
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2. The Basic Problem 


If N denotes the geold undulation with respect to the mean earth ellipsoid, 
and Ag is the geoidal gravity anomaly, then, in spherical approximation, Stokes' 
Integral reads (see Heiskanen and Moritz, 1067, pp, 92-94)$ 

N(S,A) " ■jfi^JJs(cos^)Ag(0',X’) da (1) 

R is the radius of the sphere which approximates the geold; y is a mean value of 
normal gravity on this sphere; 9,X are spherical coordinates; ^ Is the central 
angle between points (9,X) and <d',X'); and a denotes the unit sphere. Neither 
N nor Ag has zero- and first-degree harmonics under the assumption that the mass 
of the reference ellipsoid equals the earth's mass, Umt the normal potential on the 
ellipsoid equals the gravity potential on the geoid, and that the ellipsoid's center 
coincides with the earth's center of mass. The kernel S(cos0), known as Stokes' 
ftinction is ejqjandable in series form as 


S(y) « 

nara * 

(2) 

where cos (I* and its abbreviation 


y » cos 

(3) 


will be used interchangeably throughout all derivations. The coefficients 

tn = J's(y) P„(y) dy , n > 0 
or 

K == , n S 2 ; to = ti = 0 ( 4 ) 

are the Fourier coefficients of S(y) (i.e. when expanded in terms of the (orthogonal) 
Legendre polynomials P (y)). Figure 1 shows S(cos (i)) as a function of (I*. 

If gravity data are not available or sufficiently accurate over the entire globe, 
then the integration in (1) may be limited to a spherical cap a centered at the 
computation point. This procedure is associated with an error; specifically 

6Ni = J Js(cos (&)Ag da - Tfr J Js(cos f )Ag da (5) 

a ac 

a 

(Technically, this is the negative of the error, but most authors seem to prefer the 
formulation in (5).) The error kernel AKi in this case is (see Figure 2) 
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S(cos ^ ) 




aKi(cos#) 



( 7 ) 


. . fO • 0 < ^ s* lib 

4Ki(co« $) - ii> 0 < $ < n 

where is the radius (spherical distance) of the cap - 

Expanding the surface function Ag into spherical harmonic components 
» we may write briefly 

* E (8) 

ApK 

The spherical ha;?monic functions Ag„. are eigenfunctions of the integral operator 
of the type 

•) S(cos <D) da ( 9 ) 

That is, this operator applied to Agm returns Agn. multiplied by a constant, or 
eigenvalue (depending in this case oidy on the degree n); 

J J S(cos 0) ..'a « 2rttnAg.,,(8,X)} n,m5*0 (10) 

V 

This useful fact is readily established as follows. Explicitly, we have 

* Afl, R|,,(8,X) + Bjii Sni(8,X) ; n,m&0 ( 11 ) 

where 'i^i, ^ are fully normalized spherical harmonic functions, and Ani» 

Bn, are constant coefficients (see Helskanen and Moritz, 1967, p. 33). The 
series expansion of S (equation (2)) is now substituted into (10), and using the 
addition theorem for Pn(costii) (ibid, p. 33); 

Pn(co8 $) « 0 (12) 

we obtain 

JjAg„.(0'.V)S(cos^i))da « E ^J7lAn.R;:.(e’,X»)+Bn.ir-(fi'.A»)l • 

r 

• E [Rp,(«.A)Rp,(e',X‘) +Sr,(0,X) Sr,(0',X')l da ; n,m& 0 (13) 

q xso 

Invoking now the orthogonality of spherical harmonic functions, we get 

J Js(cos0)Agn,(0»,X’)da *|t„ lAn.Rn,(0,X)JjR®,(e;X)da+Bn,s;40,X)JJs„®(0’,X’)da] 
a a '‘a 

= |t„» 4nAg„,(0,X) = 2ntnAg,i,(0,X); n,m^0 (14) 

thus proving equation (10), This result holds in general for any kernel function 
depending only on $ and will be used repeatedly. It conveniently provides for 
the immediate conversion of an integral over the sphere a to a series of harmonics. 
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( 15 ) 


Therefore, i'e expansion of the error 6N« into spherical harmonics is 
simply achieved by first expanding the error kernel (equation (7)); 

<^Ki(y) - 

nmO * 

where y»cos^ and the Q^n are the Fourier coefficients of this ex})ansioni 

Qin - |.,^Ki(y) Pn(y)dy (16) 


J_’^'s(y) P„(y)dy 


(17) 


Here, yp »cos (tin. As above, we note that the constants 2irQ,n arc the eigen- 
values of the integral operator in (6), Then substituting (8) Into (6) yields 


or 


6Ni 


r A 




2y 


* -I; E Qln^gr, 

il SS3 


4S^ 


(18) 


Where 6g„ « E Ag„» 


■ kC> 


(IB) 


(Recall that by assumption Ago^* 0; also Agi -0.) The cocfCicients Qi«, being 
functions of i&o , are generally known as Molodcnsldi's truncation coefficients (or 
functions). The rigorous evaluation of Q^„ according to (17) has bc'cn undertaken, 
e.g. , by Paul (1973) who developed an accurate recursion formula (sec also 
Hagiv/ara, 1976). Hsu Houtze and Zhu Zhuowen (1979) have derived approximate 
formulas for w'hen n Is large; see also Gancko (1977). Thus, the error In 
the undulation comnritted by neglecting gravity .(inomalles Ag outside the cap ext 
can be computed according to (18). Usually, fur numerical studies such as in this 
report, the error is estimated by a global average of 6 Ni , i.e. the root mean 
square (RMS) value ; 

R^ 


(6Ni)‘ 


4r 


E Ql1.Cn 


nsS 


( 20 ) 


where the quantities c„ are degree variances of Ag, given on the geoid (sphere of 
radius R). The derivation of (20) may be found in Heiskanen and Moritsu (1967, p. 
261). The degree variances c^, not known to a very high degree, are often evalu- 
ated according to a model (see Tscherning and Rapp, 1974). 
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Formula (20) has been applied frequently in practice, for example, by Kapp 
and Rummcl (1975) who used a slightly modified version, taking into account a 
higher-degree reference field (see section 3){ also Wong and Gore (I960) and 
Fell (1978) have investigated its characteristics, in each instance, it was observed 
that the RMS error, paradoxically, does not always decrease as a larger field of 
anomalies is Integrated (cap sIko is increased) (see Figure 5). Mathematically, 
this phenomenon is easily exijlalned. Recall that is a function of , and If 
yo * cos 00 , then 


rl /I 

^ Qir. ’ - Sin lilo Ip- Qjn , n > 0, » 0 (21) 

Assuming that the (imetion Qi„ has a minimum value, it must occur when Q;„~ 0. 
From (17) 


Qin « vS( yo) Pnf yo) , n 0 , yo I (22) 

This is 550 ix> for all n if S(yo)*0. Consequently, If the RMS error 6N, (equation 
(20)) has minima, they are at the zeros of S(cos ib), namely at (bo « 38?962073 and 
ifto * 117? 66153, as seen In Figure 5. This strong dependence of the error on the 
kernel is the motivation for the discussions of Meissl (1971b) and is the basis for 
any modification to the kernel, the djjective being to reduce the error. 


In summary then, tlio geold undulation, formally the result of an integration 
over the entire sphere, is computed by limiting the Integration to a spherical cap, 
thus incurring an error. The following sections review the modifications to the 
error kernel as proposed by various authors. The error kernel is, in theory, 
modified only in its definition over the cap, requiring of course a corresponding 
change in the kernel of the integral of gravity anomalies in order to preserve the 
integrity of the original Stokes’ formula (1). 

A strong analogy exists between this idea and the concept of window functions 
in sjxjctral analysis. By restricting the integration to a cap, we ’’see" merely a 
portion of the global gravity Information, as If through a window. Thus, a change 
in the kernel, l.e. in the weights of the integration (applying a different window 
function; e.g. giving less weight to Uie Ag near the edge of the cap) can, under 
certain circumstances, lead to a better approximation of the true undulation. 


3. Melasl*g Modifloatton 


The following formulas, which form part of Meissl's (1971b) derivation of a 
recurrence formula for the coefficients, hold for any kernel function. Here, 
they are specialized to Stokes' htnction. 


Let f, g be two hinctlons, twice differentiable and square integrabie; then 
the derivation of the following analogy to Green's second identity may be found in 
Meissl (1971a, p. 43); 


(fv^g - g7®f) da = J(fvg - gvf) 'C dOB (23) 

where B is any sub-area of the unit sphere a. 6 b is the boundary of B, and 
u is a unit vector tangent to B and normal to 6 b such that it is directed away 
from B (see Figure 3). 7® and ^ are the Laplacian and gradient operators, 
respectively, both formulated in the spherical coordinate system r, 9, X* holding 
r constant (»1), For example, 


.,3, _ ?^“f . I S®f . 1 ?i®f . 2 


2 ^ cote ^ 

r* :^e 


which if r = constant = 1 becomes 


7®(f|rBi) “ 


1 a^f .aht 


(24) 


(25) 


Figured. Region B on Unit Sphere a 



The reformulation of equation (17) for Q^„ uaing Green's second identity 
as expressed in (23) will enable us to stody the asymptotic behavior of these co- 
efficients. Without loss in generality, the pole of the spherical coordinate system 
may be rotated to coincide with the point at which the undulation is computed. 

Then ^ is the co -latitude and both S(cos0) and P^cos^) are independent of 
the second variable, say the azimuth ot . In equation (23), we may identify f(i/i ,a) 
with S(cos^), v®g(0,a) with P„(cos^), B with a - O: , and 6 b with which 
denotes the circle bounding the cap a - Oh • Then (23) becomes 


f [ [S(cos(l)) Pn(cos0) -7~^Pn(cos(^) S(cos(())) Sint/) d^ da 

^ 26 ) 

= f lS(cos^)V(V”'’p„(cos0))-7-® P„(cos^)vS(coB^^))l •i)d)9 

where V”® is the inverse Laplacian operator, i.e. (V® f) ® f = (V“® f). 

Now with y ~ cos ll> , (25) is transformed into (with ^ ^6, X a a ) 



(27) 


Since S(cos0) and P„(cos(^)) are independent of a, the second term in (27) 
vanishes. The following are familiar formulas for P„(y) (Hobson, 1965, pp. 32-33): 


(n + l)P„^,(y) +nP„„i(y) = (2n+l)y P„(y), 

n > 1 

(28) 

(1-y^ R'(y) = n(P„.i(y) -y I?.(y)) 

t 

n>. 1 

(29) 

y P„'(y)-lV_i(y) = nP„(y) 

1 

n'=i 1 

(30) 

(2n + l)P„(y) = P„lx(y)-Pnli(y) 

9 

nfe 1 

(31) 

(the primes denote differentiation with respect to the argument y ). 

Using (29), 

and (27), it is easily verified that 




v®P„(y) = -n(n + l) P„(y) 

9 

n> 0 

(32) 

so that V ^p„(y) - P„(y) 

f 

n 1 

(33) 

The boundary 0 is the circle 0 = (I)q with (linear) radius 

sin^o ; and hence the 


differential arc element is d/9 = sin(i)o da , where 0 ^ a ^ 2n . Also, the unit 
vector u is in the direction of decreasing (i). Implying that Vf*C, being the 
directional derivative of f in the direction of i), is simply 
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V f • VJ « 


’i-L ^ 

— and — ^ 


» + sin f'(y) 


(34) 


btl> “““ rsil) 

With these considerations, (26) reduces to 

9 1 (*3TT 

[S(y)P,(y) P„(y)v»S(y)Jdy = ;^^JjS(yo)P„7yo)-Pn(yc^S’(yo))(l-y|)da 

(35) 


~2rr( 1-Vr?) 

■■■n(n + l) lS(yo)P„’(yo) - P„(yo)S»(yo)] , n> I 


where the operator V® is given by (27), Our original aim to reformulate the 
expression for Qi„ is achieved by recalling (17); 

“ S^i°P"(y)'=^®S(y)dy-JJj||;^tS(yo)P»^ n& 1(36) 

This, in essence, is equation (B.50) of Meissl (l&71b). 

Now if n , as Meissl notes, the first and third terms of (36) approach 
zero as 1/n^ while due to Pn'(y) such rapid convergence is not guaranteed for 
the second term (see (29)). Therefore, the condition S(yo) = 0 accelerates the 
approach of the coefficients Qj„ to zero and the convergence rate of the series 
for 6 Nx (equation (20)) is increased. This observation is the basis for Meissl’s 
strategy to reduce the error . 


From equations (1) and (5), we may write 

N = J JtS(cos|()) -SoJAgda + 6 Ni JjSoAgda 

where So - S(cos<i)o) and where the error is now 

fiNg = 6 Nx + J Jso Ag d(T = J JaKo(cos V) Ag do 
CTc V 

with the error kernel (see Figure 4) 

AKa(cos ti)) = * 0 < l/) 5.' iilo 

L S(cos 0) , «/)o < 0 S n 

AK 2 (cosiii) can be expanded in a series of Legendre polynomials; 


^Ks(y) = Z 


2n + l 


Qsn Pn(y) 


nsso 2 

where the Qg„ are the Fourier coefficients; 


Qsn = J’^AK 2 (y) P„(y) dy , n :> 0 

Y 

f [S(y) - vSoJp,,(y) dy 
•'-1 


(37) 


(38) 


(39) 


(40) 


(41) 


r X (.>h 

So Pn(y) dy + 
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aK2(cos^ ) 



n s: 1 


(42) 


Qan « 0 + rtS(y)-SolP„(y)dy, 

•'-i 

which follows by tha orthogonality of Legendre polynomials. Proceeding as before, 
the constants 2 ttQ 8„ are the eigenvalues of the integral operator in (38), and the 
RMS error ^3 is given by 

f (13) 

^ V U n sbS ^ 

Applying the kernel Ka(y) =S(y) - So to equation (35) and considering (42), we 
obtain 

== iJiSTTi) r^’P«(y)'^®S(y)dy--j;i^ „ ^ 1 ^44^ 

since 7*(S(y)-So)»V*S(y), S(yo)-So*»0, and a7lS(y) -Sol |j,„j,^=S»(yo). Therefore 
as n-»®, the coefficients Q3,, , lacking the second temi in (36), approach zero 
faster than the Qi„. A comparison of the hinctions AKi(y) and AK^ty) (Figures 
2 and 4) makes this equally evident. AKg(y) is a continuous function, while AKj(y) 
has a jump-discontinuity at y = yo. Since each of the partial sums of the expansion 
(15) of AKi is a continuous function, the series converges rather hesitantly to the 
discontinuous function. Particularly in the neighborhood of yo (if yo is not a zero 
of S(y)), the partial, sums exhibit considerable oscillations in their attempt to 
accommodate the discontinuity (Gibbs phenomenon - see any textbook on Fourier 
series) and in the limit converge to So/2 at y = yo . It is therefore also obvious 
that the series for the continuous function AKg enjoys a much Improved rate of 
convergence. This rate is further enhanced by removing also the discontinuity of 
the first (and higher) derivative at y =* yoj however, this is not pursued here (see 
Meissl, 1971b, p. 52; cf, Hsu Houtze and Zhu Zhuowen, 1979). It must be noted 
that the convergence rate is not a critical issue in the computation of the RMS error 
with today’s computers. The upper limit of the sum (equation (20) or (43)) is easily 
taken to 2000 or 3000 which Is normally sufficient for at least millimeter accairacy 
(somewhat less for Vo =0°) in the approximation of the error series by a finite sum. 

Although the coefficients Qg„ converge to zero more rapidly than Qi„, this 
does not ensure an RMS error that is smaller than 3^(see equations (20) and 
(43)). Infect, for cap radii less than 40®, the error 6Ng exceeds 6 Ni , as 
clearly seen in Figure 5. Table 1 lists values of Qsn versus Qi„ for 4>o = 10® 
verifying the accelerated convergence of Qg^ to zero, but also showing that more 
’’power” is shifted to the low-degree harmonics of the error kernel AKg » whence 
the larger error. Evidently, the error SNs is significantly smaller than 6 Ni for 
a value of $0 such as 65®; but considering the current state of world-wide gravity 
data, caps with radii larger than 20® or 30® are usually deficient in coverage that is 
both dense and accurate. 
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30.0 4 Figure 5. RMS error in N due to truncation at cap ladhis it>o 

using die usual Stokes’ fiinction (imrve 1) and 

I \\ Meissl’s modification (curve II). 



ap Radius 


Table 1. Truncation CoefCicienta of the Error Kernels 
AKi, AKa.AKa #o » W®) 


n 

Qln 

Qan 

Q3„(m = 20) 

0 

-.414 

-.201 

.034 

1 

-.411 


-.051 

2 

1.593 

1.801 

.032 

3 

.599 

.802 

.030 

4 

.274 

.471 

.027 

5 

.118 

.307 

.025 

6 

.030 

.210 

.021 

7 

-.023 

.147 

.017 

8 

-.056 

.103 

.013 

9 

-.076 

.072 

8. 36 X 10"® 

10 

-.086 

.049 

3.46 X 10“® 

15 

-.073 

-3.65 X 10"® 

-.023 

20 

-.025 

-.012 

-.047 

25 

.013 

-7. 80 X 10"® 

.020 

30 

.026 

-1. 65 X 10"® 

4.38 X lOr* 

50 

-.013 

-1. 33 X 10"* 

2.05 xior® 

100 

3.89 X lO"® 

-1. 54 X lO"* 

-3. 79 X 10'* 

ISO 

-7.f 6 X 10"^" 

9. 89 X 10"® 

-1. 78 X 10“® 

200 

-5.91 X 10'* 

-4. 73 X 10"® 

1.41 X 10"* 

300 

-8.77 X 10"* 

3. 51 X 10”® 

1.22 X 10”* 

500 

4. 10 X 10"* 

8. 48 X lO"*^ 

-6.00 X 10”® 

1000 

1.27 X 10"* 

-4. 61 X 10"^ 

-1,78 X 10"® 

1500 

2.72 X 10"® 

-3. 13 X lO""^ 

-3.55 X 10"® 


Meissl's modification of the kernel therefore has practical applicability if the iow- 
degree harmonic components of the gravity field are known (not necessarily perfectly); 
the corresponding terms then do not contribute to the truncation error. Let m denote 
the maximum degree of the available potential harmonic coefficients. Then one may 
compute according to the conventional approach 


A 

Ni 

with an error 


R 

4n y 


J J S(cos 0) Ag da+ A Ag„ 


(45) 


^ E Qin e(Agn) + f- 

(assuming no errors in the data of the cap Oc 
the truncated kernel: 


E Qin Agn (46) 


= ■+1 

), or by removing the discontinuity of 


- 15 - 








(41) 


A 

Na 


With an error 

6Na' 


9 

JJlS(co8 ii»)-Sol Agda +*^ Eg Qa" 


* OD 

^ E Qan ^ E Qan ^fi>n 

nPtS ms bH 


(48) 


The quantities Ag„ are the harmonic functions computed from the coefficients. The 
corresponding error is 


6(Ag„) = Ag„ - A g„ , 2 s n s: m (49) 

The RMS estimates of 6 Ni' and ONg are given by 

6Nx' = E Q^„® 6c„ + E Qin® c„f (50) 

l-n**3 nsi«+i "J 

and = ^FE Qan® &« + E (51) 

n=*+l ■ 


where dc„ is the degree variance of the errors in the harmonic components of 
degree n. The validity of (SO) and (51) follows immediately under the hypothesis 
of zero correlation between fte errors in the coefficients and the harmonics above 
degree m. 


Wong and Gore (1969) adopted a conceptually different method (with no claims 
to reduce the error) in which the first m - 1 harmonic components are subtracted 
from the kernel (see also Fell (1978)). Recalling (2), let 

, , r> 2n + l 

S (cos 0) = 2 j 


2 t„P„(COS0) 
nsii+1 ^ 

then the computation of the undulation according to 

^ " 4^ JJs'(«os0)Agdcr 

is associated with an error 

6N3 = ‘T^r JJs(cos0) Ag d(7 - JJs*(cos0) Ag dcrj 

^ (7 Oe 

= 4f^[JJs(cos0)^Eg^g„da+ ^Js(cos0) Ag- da + 

+ J Js*(cos0) Agda - JJs"(cos0) Ag da] 
a-ofc ^ 

~ 4^[ J J Js(cos0) J Ag„da + 

a ^ a 


(52) 


(53) 


(54) 


JJ[S(cos0)-Sl(cos0)lAg"da - JJs"(cos0) E^Ag^da] 
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where 


(55) 


Ag* ■ E ^iPB 

The laet two terms in (54) vanish by the orthogonality of harmonic ftinotions, since 
S(cos||)) - S”(cos(t>) possesses no harmonics beyond degree m and S* has none 
below degree m + 1 . Equation (54) can therefore be rewritten as 


+ J* Js^cos|i»)^Eg^g„da+ JJs(cosiii) E^Ag„do] (56) 


Given potential coefficients to degree m (subject to errors ^(Ag,^)), the second and 
third terms of (56) may be included in the evaluation of the undulation. Hence we 
calculate 


Na “ 4^1 jS'(cos(|i)AgdO+ ^ E (Q 3 „+t„) a|c, 
<7c 


With an error 


^ E (Qa. + t„)6(Ag„) Qg„ Ag„ 

“An b3 “ 7 


Where 


(57) 


(58) 


(59) 


Qsn = f S-(y)P„(y)dy , n> 0 
are the Fourier coefficients of die expansion of the error kernel (cf. equation (56)); 

AK3(cos 0) = I ® / V t ^ (60) 

' lS*(cosi&), ' 


and the coefficients t„ are given by (4). Numerical values of Q 3 „ for it>o = 10° 
are listed in Table 1. The RMS error Ssjt from (58), is 


6N. 


E (Q» + 6cn + E Qr/ 0 J 

“7i-ns2 r. p:«tX J 


(61) 


Note that AKs » expanded in the interval [-1,1] as a series of Legendre polynomials, 
contains all harmonics from degree zero to infinity. The constants 2 ttQ 3„ are also 
the eigenvalues of the integral operator in (56), whence the results (57) and (58) 
follow easily (see section 2). 

This section is concluded with a numerical comparison of the three alternative 
RMS errors in the undulation that are caused by the lack of complete and global gravity 
data, as well as the errors in the harmonic coefficients. Specifically, these error 
estimates are given by equations (50), (51), and (61). Equation (50) is the RMS 
error if Stokes' function S(cos((>) is applied to the cap Oc » equation (51) corres- 
ponds to the use of S(cos 0) - So in Ob ; while (61) represents the error when 
S *(cos((>) , as defined in (52), is used. 
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The gravity inomaly degree variances are modelled according to Tscherning 
and Rapp (1974): 


ovf n a- 3 

" (n>2)(n-f24) ^ 


(62) 


where s ■ .999617; and the degece variances 6o„ are determined from 

flc„ - y»(n-l)= , ns 2 (63) 

where the are the degree variances of the errors fiCm t ^r« in hilly normalized 
harmonic coefficients of degree n . The formula for 4n 


4« - y; C(fiCr.)*‘i-(^S„)3|. n s 2 

■TiO 


(64) 


The estimates of , as listed in Table 2^ were obtained by substituting into (64) 
the standard deviations of the coefficients of the GEM 9 (20,20) solution (Lerch etal. , 
1977), 


Table 2. Estimates of Degree Variances of GEM 9 
Potential Coefficient Errors 


n 

<n X 10'® 

n 

L X iC® 

2 

.000037 

12 

.004471 

3 

. 000341 

13 

. 006040 

4 

. 000205 

14 

. 005048 

5 

. 000964 

15 

.006343 

6 

.000656 

16 

.005538 

7 

.002130 

17 

.006873 

8 

.001522 

18 

.007575 

9 

. 003536 

19 

.006882 

10 

11 

. 002940 
. 005949 

20 

.007038 


The Molodenskii truncation coefficients can be computed, for example, 
by the recursive formulas of Paul (1973). Then with relatively little additional 
effort the coefficients Qsn are obtained from (36), (44), and (29) as 

Qsn = Qxn +-j^[Pn-i(yo) -yoP„(yo)l. n S 1 (65) 

where So = S(cos0g), yo = cos 0o. Finally, the coefficients can also be 
modified to yield the coefficients Qa„. Substituting (52) into (59) and considering 
(2), we have 
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n 0 


( 66 ) 


Q» • Qi» - E Ij(y) P«(y) dy. 

rma « *'-1 

- Q» - £ e,. 

r*a ^ 

where 

e„ - f’^°Pr(y) P«(y) dy, r.nSO (67) 

Then Hobson (1965i p« 36) gives 

Or„ - ^V.|j(rTn4l) ^-i(yo) - r P„(y«^ Pr- i(yo) +yo(r-n)P«(yo) Pr(Vo)l, (68) 

Also , r^n. r,n>0 

%n" J^lPr^i(ye)-IVM(yo)l. r>0 (69) 

©0.0 * 1 + yo 

Or.r - 5^((2r-l)er-i,r-l+yo(lf(yo) + Pril(yo))-2Pr(yo)Pr-l(yo)l. r>0 (70) 

Equation (09) is obtained by inserting (31) into (67) and noting that Pr(-iV » (-l)*^. 
Equation (70) is derived in Appendix A. The values R * 6371 km and y » kM/R® * 
982026. mgal (where kM 398601 km®/sec® is the product of the gravitational 
constant and the earth's mass) were used in all computations. 


Figure 6 shows the RMS errors for cap radii 0o “ 0° , . , , , 40® with the assump- 
tion of no errors in the harmonic coefficients up to degree m » 20; i.e. , 6c, , » 0, 
n " 2, ... ,20. Clearly exhibited is the substantial difference in magnitude between 
the SNi' and SNs' truncation errors, even for smaller caps. For example, at 
>l>o= 10®, - .82 m, 8Na' - >26 m, and » .82 m. Doth the unmodified 

error kernel and the kernel of Wong and Gore (1969) produce similar truncation 
errors; but the introduction of errors into the reference field causes the two 
corresponding RMS errors to diverge considerably, as seen in Figure 7. This 
characteristic Is attributable to the coefficients t;, in equation (61) which are 
independent of the cap radius ifio ; and therefore, Sn,' reflects primarily the 
large error due to the erroneous harmonic coefficients (e.g. , an RMS contribution 
of 1.61 m for (l)o “ 20®). 

Of the three alternative methods, considered in this section, to compute un- 
dulations, Figure 7 evidently favors (for i(io ^ 4®) the simple modification elaborated 
upon by Melssl (1971b). Furthermore, the feasibility of Integrating caps of radius 
between 5® and 10®, and yet achieving half-meter accuracy (assuming a 20-degree 
reference field such as GEM 9 and neglecting integration and gravity data errors), 
is manifestly demonstrated. Figure 7 gives, for example at 0o = 10°, 5 Ni = 1.09 m, 
llNg = 0. 41 m, and iSlg = 1,67 m. It may be noted, again, that the Improved 
convergence rate of does not produce a corresponding reduced magnitude for 

^0 < 4 °. 
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RMS Error (meters) 



,o 


21 


lO 



Cap Radius ^ 





Gap Radius 





4. Molod0nBkli*B Improvement 


Molodenskii's a|)proach to reduce the truncation error, originally (Molodenskii, 
1958) devised as the determination of the most rapidly converging error series, adds 
sophistication to the concept of section 3 as more than a mere constant term Is sub- 
tracted from Stokes* function. Recapitulating the problem, the error in the undulation 
derived from an iutogration of gravity anomalies over a spherical cap (radius t'o)# 
using Stokes' function as the kernel, may be represented by an integral over the entire 
sphere where the error kernel is zero over the cap (see section 2), When this error 
Is expanded in a series, the coefficients converge to zero slowly as a result of the 
l.emel'B discontinuity at lb (bo . However, if the error kernel were to be redesigned 
so as to eliminate not oiUy its discontinuity, but also the discontinuities at ((lo of all 
its derivatives, then we should have the smoothest error kernel possible with the 
most rapidly converging error series. This modification cf the error kernel amounts 
to an analytic continuation over the cap of the truncated Stakes’ function (l.e. , from 
0*0 to 0 « But then we would encounter a dilemma since Stokes* function Itself 
is analytic for 0 < (I) ^ n, and the error kernel continued thus would coincide with 
Stokes' function exactly for 0 0 sJ n . That is, an analytic function is determined 

uniquely by the values of it and all its derivatives at a single point - in our ease, these 
values are already specified for all points in the interval (0o i ”I • Consequently, as 
easily recognized from equations (71) and (72) below, the kernel for the integration 
over the cap would be zero, nullifying the contribution of the gravity anomalies in Ob . 

A viable alternative (and this is how Molodenskli proceeded) is to approximate 
Stokes' function in the Interval (0o i "1 by a finite sum of Legendre polynomials. 

This new kernel enjoys all the properties of being analytic over the entire Interval 
(0, tt] , but an additional error is introduced since it does not coincide exactly with 
Stoksis' function for 0 q 1 0 < tt . Specifically, let the undulation be computed according 

to A R r r 

N « ITTy J J Km (COS0) Ag do (71) 

Oc 

where Km(cos0) = S(cos0) -'S‘h(cos0), 0 < 0 ^ 0o (72) 

and the function 

Sii(y) = E ^^|^s„P„(y) (73) 

is to be defined by the coefficients so as to render the "best" approximation to 
S(y) for 00 ^ 0 ^ rr (see Figure 8). This is not to be confused with the approach 
followed by Wong and Gore (1969) since s„ tn - the concern here is only^with the 
interval [0o,tt] (see below for the determination of s„). The error of N is 

N-A = J Js(cos0)Agda - jj^||lS(cosi/)) -§'s-(cos0))Ag da 

' O CTc 

-A-f '’s(cos0)Agda + rr§'ir(cos0)Ag da (74) 

4nyJ^yJ 
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OF POOR QUAI^ra 


Fig*4re 8. Stokes' ftmctloD vs. its sn>roximstion 
3V t n ■ 5 » $0 • 20 *^, 
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A 

N - N 


B 


J[ S(C 08 di) - §V(cos v*)! Ag dff 

'O - Oj 



(76) 


The last equality follows by the orthogonality of harmonic furetions, since Sji(y) 
contains no harmonics above degree U . The fir^t term in (75) represents the 
error due to the discrepancy l^tween S(y) and Sir(y) in rrj . (If n ", 
then this term is zero since S|r(y) -♦ S(y>; but this is true for all y in [-1,1), 
So that also N 0.) If potential coefficients to degree m are available (with 
possible errors), then we may set n « m , thereby essentially eliminating the 
second term of the error (equation (75)). Thus, the undulation is computed accor 
ding to 


A 

Nm 

with an error 


= R 

4Try 


f J { S(cos \li) - §',(C0B (|>» Ag da + ~ j; s„ Ag„ (76) 

Oc = 3 


-S,(cos(&))Agda + ^ f; s„ 6(Ag„) (77) 

where, obviously by (73), the constants 2nSn are also the eigenvalues of the second 
integral operator in (75). 


Before^proceeding to ftirther simplify (77), S, must be determined explicitly. 
Recall that being a finite sum of polynomials, is to approximate S as closely 
as possible in the interval l(&o» • le the recent papers which explore this method 

of reducing the error (Dickson, 1979; Fell and Karaska, 1979; cf. Hsu Houtze and 
Zhu Zhuowen, 1979) the function S, is determined through a "least squares adjust- 
ment". More precisely, by quantifying the concept of closeness through the definition 
of a norm, a well known and powerfiil theorem of Fourier analysis yields the desired 
result. Let f(x) be a function, square integrable in the interval [-1,1); then the 
norm of f may be defined as the square root of 

llfp = J_'^f"(x)dx (78) 

Now change the variable y (= cos (li) to 


y = kx + k-1; 


k = cos'* 

2 


h ( 1 + cos;l*o) 


(79) 


As y ranges from -1 to yo, x varies from -1 to 1. We note that while the 
Legendre polynomials Pn(y) do not form an orthogonal basis in the interval (-1, yo), 
such a basis is formed by the polynomials Pr(x) . ^,(y ) given by (73) (with ti = m) 
is a polynomial of degree m and can equally well be expanded in terms of P„( x): 


S.(y) 



2r-H 

2 


Ur Pr(X) 


(80) 


where, formally, the Fourier coefficients u, are given by 
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Ur - J^^S,<kx+lc-l)P,(x) dx - S,(y)P,(^~^)dy} k^O, rfeO ( 81 ) 

The *’beet” approximation of Sg to S In the interval [-1, yol» or equivalently for 
-1 5x^1, is then obtained by minimizing the norm of the difference: 

J|(S(y) -S^,(y)l®dx - min, 

or r*tS(y) -S,(y)l®dy - min. 

**-i 

The minimum condition (82) is fulfilled if (Davis, 1975, Theorem 8.51, p. 171) 
the coefficients of Sa(kx+k- 1) (a finite sum of polynomials Pr( x) ) are the Fourier 
coefficients of S(kx-fk-l); 

Ur •- r S(y)Pr(x)dx (84) 


(82) 

(83) 


Considering equation (77), let 


4K,(cos*;, = {s(e„s*)-s.(co.*): 

which may bo expanded in a series 

y 2n + l 

AKM(y) = L -T“<^r.Pn(y) 

n “ 0 * 

The coefficients are given by 

Qmh = j °(S(y) - S.(y)l P„(y) dy 
•J-i 

= Qln “ Qn 


0 S l(fo 


(85) 


( 86 ) 


(S7) 

( 88 ) 


Qn = J"°S.(y) P„(y) dy (89) 

and the are the usual Molodenskii truncation coefficients given by (17). In 
view of (84) and (80), §^a(y) is the truncated expansion of S(y) in the interval 
[-1» yol • The immediate consequence is that 

Qnn = 0 , 0 s n iS m (90) 

Indeed, from (87), (80), and (79), and for n ^ m, we have 

Qin-Qn “hf [S(kx+k-l) - S.(kx+k-l)l Pn(kx+k-l) dx 

CO 

= k fV E Ur Pr(x) ) P„(kx+k-l) dx (91) 
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Now P„(kx+k-'l) is a polynomial in x of degree n s m, which is also expressible 
as a finite sum of I.egendre polynomials P*<x), 0 s: s n (see section 6). Therefore 
P„( kx +k - 1) contains no harmonics above degree m ; also S( kx +k - 1) - S,( kx +x - 1) 
has none below degree m (i.e, , r s m + 1 in (91)) - the orthogonality of Legendre 
polynomicals ensures the result (90). The remaining coefficients (n > m), 
characterizing the error (see below), are expected to be relatively small in magni- 
tude, this depending also on the cap radius $q. 


The error 6 Nm' (equation (77)) is now more suitably written as 

■ 00 

«N.’ = ^ E ^ E Q.„ Ag. 

^7 n=a *7 n = ii+l 

The corresponding RMS error is 

®Nh' = E 6c„ E 

-57 n=s3 


n = *+l 


J 


(92) 


(93) 


The second term in (93) arises from the difference between S^(y) and S(y) for 
-X s y s: yg. However, we cannot choose n arbitrarily large, anticipating a re- 
duction of the truncation error beyond significance, unless we have an equally large 
arsenal of known harmonic coefficients. For example, if we let n > m in equation 
(75), then without difficulty, one can deduce 

■ T 


n = 1 + 1 


nasn + l 


where Nmi is given by (76) with S„(cos0) replaced by SB(cosiii). Alternately, 
if IT < m , then it follows just as easily that 

T « 

® 4^ J 1^ ^ E^s„ Ag„ + ^ ^ ^ Be (95) 


With an RMS error 

TJ > ~ ^ 

^3' = -5^1 E 6c„+ E + E ^'nCnT' 

'^7>-n = 3 n=“+l n = «+i 

The results of a numerical comparison of (93), (94), and (96) are presented in 
section 5. 


(96) 


Assume, for the moment, that 6( Agn) = 0 , n = 2, . . . , m , so that the appli- 
cation of Schwartz's inequality (Davis, 1975, p. 134) to (77) yields 

1®^"' ^ ^(£°[S(y) -§'.(y)f dy)^(JJ(Ag)®da)^ 

^ ^ r(J_JS(y)-S.(y)l^dy) (97) 

where, for all 9,X, J J( Ag ( 0' , X') )^ dcr ^ T® = constant, and where I*! 

a-ac 
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I denotes the absolutej^alue. We see that implicit in the attempt to minimize the 

I dilference between SH(y) and S(y) in the interval [-1, yd is the minimization, 

^ in some sense, of the error 6 Nm' , The original objective to obtain the most 

rapidly converging error series is inconsequential as far as equation (93) is con- 
j cemed (If 6c„ s? o ), since the error is not characterized solely by the coefficients 

of the analytic function S]i(y). Only in (94), where n>m, do the s„ play a 
part in the truncat ion error. Nevertheless, the ix>int must be stressed that an 
increased convergence rate does not imply a decrease in the error (see section 5). 
Molodenskii et al, (1962) adopt the approach which starts v/ith equation (77) subject 
to the condition (83), and no mention of convergence rates needs to enter the dis- 
cussion. 

Interestingly, the error kernel AKm , given by equation (85) is discontinuous 
at (fi = «po ! The discontinuity is relatively "small”; nevertheless, Meissl's 
strategy is Immediately called to mind - subtract from the error kernel its value 
at (!)o and thus Increase the rate of convergence of the error series. Table 3 
shows the subsequent increase in the RMS truncation error ( 6c„ = 0, n = m “ 20), 
as the improved convergence rate is not realized until n is very large (n > 170, 
for = 10° , see Table 3a). In a differed approach, Hsu Houtze and Zhu 
Zhuowen (1979) propose that the function S, be determined by simultaneously 
removing the discontinuities of AKm and its first derivative at 0 ~ This they 
accomplish by subjecting the minimum c^dltion (83), through the use of Lagrange 
multipliers, to the constraints S(yo) » S,(yo) and S’(yo) “ s';(yo). Whether 
this method finally produces smaller truncation errors is not clear since the 
difference between S(y) and S,(y) in the interval [-1, yo)* as measured by 
the norm (78), must necessarily increase under the above constraints; and it 
was shown already (equation (75)) that the truncation error is directly influenced 
by this difference in the kinctions. 
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Table 3. RMS errore in N t aaauming errorlesa potential coefficients 
to degree m » 20, using Molodenskii's kernel and his kernel 
modified by removing its discontinuity^. (Errors in meters.) 


^0 

W 


1® 

1.93 

2.53 

2® 

1.13 

1.74 

6® 

.28 

.47 

10® 

.03 

.05 


* using the error kernel AKwaiy) “ {s(yf ' l^jy*)^! ^ J y 


Table 3a. Coefficients and Qhs of the error kernels 
AKm and AKm 3 (see Table 3)j “ 10® . 


n 

Qm 

Qm,3 

30 

100 

200 

300 

1500 

-5.98 X 10"* 
4.97 X 10”® 
-2. 50 X 10"® 
-2.05 X 10"® 
6.00 X 10” 

-1. 28 X 10-® 
-5. 37 X 10"® 
-1. 15 X 10"® 
1.29 X 10"® 
-8.11 X 10"® 


< 1 

< yo 




5, Computational Procedures and Reauits 


The previous section dealt primarily with Molodenskii's theory on the reduction 
of the truncation error when the undulation is obtained by integrating Ag over a cap 
Oc end from harmonic coefficients, bi this section, the theory is implemented and 
all necessary working formulas are derived. The essential quantities to be deter- 
mined are the coefficiei^ts s„ of the function Sir(y), since they are required in 
both the calculation of Nm (equation (76)) and the estimation of ^e error (equation 
(93)). The e;q)ansions for §r(y) nre repeated for convenience: 

SK(y) = E P«(y) ( 98 ) 

nmO « 


S-(y) = t -^^UrPr(X) 

r« 0 “ 


(99) 

^here v - k + 1 

y=kx + k-l, X = , k 

cos®|“ 

(100) 

The Fourier coefficients Qm, Q,, s„, Ur are given by 




(101) 

Qr = J jSir(y) P„(y) dy , n s* 0 


(102) 

Sn = J_^§r(y) Pn(y) tly » O^n^Jn 


(103) 

Ur = J ^S(y) Pr(x) dx , 0 £ r ^ n 


(104) 

Changing variables from x to y in (104) yields 




n) , 0 r ^ n 

(105) 


Now, Pp((y-k + l)A) is a polynomial In y of degree r. Furthermore, it is 
defined for all y, in particular in the interval (-1,11, Therefore, Pr((y -k+ l)/k) 
can be expanded ao a linear combination of the orthogonal (independent) polynomials 
Pn(y) » 0 s ns r , which generate the space of all (real) polynomials of degree r 
defined in [ -1, IJ; 




hrn P„(y) 


(106) 


where the coefficients hm are given by 
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• r,n i 0 

We note that by the orthogonality of Legendre polynomials, 


h,„ ■ 0 if n > r 


( 108 ) 


since P,((y-k4l)/k) has no components beyond degree r. Substituting (106) into 
(lOS) results in 

Ur « £ ~^h,„pS(y)P„(y)dy 

K n«0 ^ •'-1 


. 1 I l^h„ Q„ . 0 s r < n 

Using the expansion (99) in equation (103), we obtain 

IT 

V 2 r-f 1 . ft ^ - 

“ w II j. hpjv ^ 0 ^ w ^ ^ 

rs n " 

where (108) has been used. Finally, from (98) and (102) 

~ r^o S 2 r+ 1 

Qn == J , L — Pn(y)dy 

r sO « 


(109) 


( 110 ) 


2r+l 


n 0 


where suitable working formulas for are given by (68) to (70). 


(Ill) 


The truncation coefficients Qi„, required in the computation of s„ and Qm„, 
have been studied extensively, and efficient and accurate recursive formulas exist in 
the literature. Hagiwara's (1976) elegant solution is used here. A recursive re- 
lationship among the hm is established using equation (31) and several integrations 
by parts (cf. Appendix B). The final result is 

k u I 2 r ■f 1 Ph r ,n— X hp,n+iT ^ ^ .. \ 

*»r+i,„ = ^ 0-^n^r+l 


■''.o' [P'*' ( k ) - P'-’ ( k )] ’ ' ^ " 

ho,o = 2; hi^o “ ' ***•' ~ "slT 

hp „ = 0 , n > r 


>(U2) 
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Formula (112) la numerically unatable (at least for ^ 30*') and is useless in 
practice. Instead, we resort to a computationally more burdensome and time 
consuming closed expression for Ivn that was developed by Molodenskii (1958): 

h,„ - (211 + 1) k’- • 0 s n < r 

p « r-n- 1, q ■ r+n 

^rr “ 2 r + l ^ * r «: 0 

where, e.g., (^)denotes the binomial coefficient Tl'/p -1)1 * Some care must 
be exercised when calculating these coefficients on a digitol computer. The com- 
putational accuracy of (113) can be tested using equation (106). This has been 
done by the author to degree r » 80 with completely satisfactory results (for 
00 - 10^, 16-dlglt accuracy at degree 0 deteriorate to 11-digit accuracy at degree 
80), 


In Figures 9 and 10 the RMS errors corresponding to Molodenskii's modifica- 
tions to the error keniel are compared only with Meissl's case. Figures 6 and 7 
already indicate strikingly that through a proper modification of the kernel the 
truncation error can be reduced significantly. Therefore, the question is whether 
to choose Meissl's almost trivial modification or to adopt the computationally more 
elaborate method of Molodenskii. Of all the alternatives shown in Figure 9 (based 
on errorless harmonic coefficients to degree 20 (GEM 9)), Molodenskii’s modifi- 
cation (i.e. n =: m ~ 20) yields the smallest truncation error for most 0o • The 
modifications defined by the condition n ^ m , as suggested in section 4 (equations 
(94) and (96), herewith n » 25 and n = 10, respectively; and 6c„=0), although 
giving truncation errors silently smaller than show no Improvement over 

the case when n_f m “ 20. From Fig ure 9 ( 6cn = 0 ) and for 0o = 10° , 6 Nh’ - . 03 m, 
6 Nhi = *09 m, 6 Nm 2 = • 15 m, and fiNg = .26 m. Introducing the errors (standard 
deviations) of the GEM 9 20-degree reference field (see section 3), but still neglecting 
integration and gravity data errors, the simple modification of Meissl for ^ 6° 
emerges as an effective means to reduce the total RMS error, surpassing even 
Molodenskii's method (see Figure 10). However, the alternative modification rep- 
resented by equations (95) and (96), with H = 10 , yields consistently the smallest 
errors for all i/>o . J'pr example, at 0o = 10°, CSms = .33 m, 3 Nb “ .41 m, 

6vl!i = . 46 m, and = . 54 m. 
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6, Numerical TeiUi with Actual Gravity Data, und the Atmomherio correction 


The computations of the geoid undulation according to the familiar un- 
modified case (equation (45)), as well as with Melssl’s modification (equation (47)) 
have been carried out with actual gravity data by Rapp (1980) for two oceanic 
areas, each spanning 30^ in both latitude and longitude; Tonga Trench Area and 
Minn Ocean Area. The (point) undulations were determined on a l^’-grid from 
an integration of l^’x I** mean gravity anomalies within a cap of radius ((lo « 10^ 

(these were determined from a combination of terrestrial and altimetric anomalies) 
and from the GEM 9 potential coefficients. Subsequently, these computed undulations 
were compared to the GEOS-3 altimeter geoid for which Rapp (1980) gives average 
accuracies of dkO.61 m (Tonga Trench Area) and i:0.71 m (hidian Ocean Area), 
Table 4 lists the corresponding mean and RMS differences between the computed 
and altimeter geoids (ibid, p. 0 and p. 9). With the same data as described above, 
Molodenskii's modification (equation (76)), as well as the alternative method with 
IT» 10 (equation (95)) were tested similarly; the results are also shown in Table 4. 
Any of the modified kernels yields a geoid that agrees better with the altimeter 
geoid by almost 50% than the geoid determined using the unaltered Stokes' function. 
The RMS differences between the geoids based on the modified kernels and the 
altimeter geoid approach the quoted accuracy of the latter. It is therefore not 
possible to rate one modified kernel better than the other from these tests. 

The set of gravity anomalies used for the integration, as mentioned above, 
was derived from both terrestrial anomalies (obtained on the actual geoid) and 
altimetric anomalies (acquired through collocation from altimetry). In determining 
the precise relationship between these two types of anomalies, one must heed the 
effect of the atmosphere (see also Rapp, 1979). The details of this relationship are 
explained below (cf. Rapp and Rummel, 1975). While these elaborations constitute 
a digression from the main text, the final results do pertain to the application of 
truncation theory. 

To simplify the discussion, the sea surface is assumed to be a stationary, 
equipotential surface, which then serves as the geoid, by definition. Moreover, we 
suppose that the geoid encloses all terrestrial masses, but not the surrounding 
atmosphere. Accordingly, the gravity potential on the geoid, Wo» is that potential 
which would actually be observed. The altimeter measurement, being purely geo- 
metric, provides directly the geoidal surface as defined above. In this respect, it 
is necessary to define only the size, shape, and position of the reference ellipsoid. 
Let the ellipsoid be centered at the geoid's center of mass and aligned with the 
earth's rotational axis. For subsequent convenience, however, the scale of the 
ellipsoid will be obtained dynamically be equating the ellipsoidal mass with the 
geoidal mass plus the mass of the atmosphere and by further defining the normal 
potential on the ellipsoid, U©, to equal W©. Given the ellipsoid's flattening and 
rotational rate, its semimajor axis is then uniquely determined (see Heiskanen and 
Moritz, 1967, p. 110). 
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Before the unduUitlons from altimetry can be proceaaed forther Info gravity 
aaomalleB (via the Inverae Stokea* formtda, or collocation - In eaaence, aohitlona 
of the boundary-value problem), theory demanda that the atmoaphere which envelope 
the geold bo removed. According to Moritz (1974), thia la conveniently achieved by 
redlatrlwting all atmoapherlc maaaos (aaaumed to be apherically aymmetrlc) above 
the elllpaold (Including fictitloua atmoapherlc maeaea between the geold and ellipsoid) 
Inside die geold. In a manner which leaves the center of mass unchanged. The intro- 
duction of mass Info the geold results In a cogeold (see Helskanen and Moritz, 1967, 
p. 141): 


N® - N-AN 


(114) 


where AN is the indirect effect produced by the change in potential 6 Wa due to this 
redistribution; 


AN 


9 «Wa 


(115) 


(116) 


From (Moritz, 1974, p. 13) we have 
. -k 

« J . (r*) 

rmt^ ' ' 

where k Is the gravitational constant, and M( r') is the mass of the atmosphere above 
the sphere of radius r*. 6 Wa is the difference between the potential due to the atmos- 
phere lying above tlie geoid and the potential of the atmosphere redistributed below the 
geold. This difference is rather small, resulting in an indirect effect, AN , of approx- 
imately 6 mm (Rummel and Rapp, 1976). The corresponding change in gravity on the 
geoid is 




rir 


(OWA)lr, 




(U7) 


that Is, the (negative) attraction of the atmosphere (which is now inside the gi'oid). 
Note that (under the simplification of spherical symmetry) the atmosphere, when 
situated above the geoid, has no effect on the geoidal gravity; and hence, tlie re- 
distribution affects gravity more substantially ( 5gA 37 mgal at the geoid, see 
Moritz (1974)). Moritz has chosen the signs such that the potential and gravity 
anomaly on the original geoidal surface become, respectively. 



W® = W.) - 6Wa 

(118) 

and 

Agf - Ags - 6gA 

(119) 


where Agg is the geoidal gravity anomaly prior to redistribution of the external 
masses. Observe that the gravity potential on the cogeoid is Wo and that the 
reference ellipsoid has undergone no changes in the process of redistributing the 
atmospheric masses. From our definitions, the normal potential on the ellipsoid 
equals the gravity potential on the cogeoid and both surfaces enclose the same mass. 
Therefore, the cogeoid undulations have no bias (Helskanen and Moritz, 1967, p. 101). 
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The Indirect effect, AN (~6 mm), In practice !• evidently neglected. The 
altimeter unduUtiona may thua be Identified as undulations of the cogeoid (with no 
external masses), having a global average equal to zero, and which therefore are 
diikcctly usable in collocation to estimate gravity anomalies. However, the resulting 
anomalies must be correctly intcrprated as anomalies on the cogeoid referring to 
an ellipsoid containing both the terrestrial and atmospheric masses. Theso altimetric 
anomalies, Ag^tt, are consistent with Ag“ (neglecting the indirect effect of (1 mm on 
the value of gravity). The terrestrial anomalies, Agr«r i on the other hand, are 
obtained on the actual geoid (with no internal atmospheric masses), but also refer 
to an ellipsoid containing the atmosphere (assume it is the reference ellipsoid defined 
above). The two types of anomalies Agr.r end Ag/i{», are rendered compatible, accor- 
ding to (119), by either subtracting the atmos|)heric correction Ag^ from Ag*^^ : 

Ag® « Agr„ - 6gA consistent with Ag^m (120) 

or adding Ag^ to the altimetric anomalies: 

Aga » Ag^n + Ag* - consistent with Agrw (121) 


At sea level, the correction Ag, is a constant, Hence, for convenience, the 
set of gravity anomalies which enter the truncated Stokes' integral was obtained by 
merging anomalies Ago (equation (121)) and Agr.r > and subsequently corrected by 
- Ag* so that the integral is validly applied. Regarding the integration, this constant 
part can be treated separately. 

Substituting (119) into any of the equations for the undulation (45), (47), (57), 
(76). or (95), we have 

" 4ny JJ K(cos 0) Ago da + A' - Ag* JJk(cos (&)da (122) 

^ Oc 

A A 

where N* represents the contribution to N from the harmonic coefficients, which 
are already based on a geoid (cogeoid) that contains the mass of the atmosphere. Let 
the last term in (122) be denoted by ANa . This is the correction to be applied to the 
computed integral (the first term). If the kernel K is the usual Stokes' function S , 
then 

S(y)dy = S(y )dy - J^°S(y)dy J = ^6g, (0-Qi,o) 

-1 -1 

= ^6gAQi,o (123) 

SimUarly, for Ka(y) == S(y)-So, using (41) 

ANa A = 2 ^ AgA j^J^(S(y) “So)dy - (Qa^o ■■ 2Sc^^ = ^Ag^^Qa^o (124) 
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And for Kn<y) *S(y)-S-(y), we obtain with (90) 

[ J (S(y) -Sr(y)) dy -Qm,o] “ ^ (0 - So- 0) 

“ *0 (125) 

For i/b “ ®Nia »» 1. 17 m, * 0.57 m, 6n«a If.ao ” 0* 43 m, and 
fiNH*lTr»io" 0,60 m. These values should be added to (R/4rf y)J^/K(cos 0)AgGdcr + N’ 
(equation (1^2)) to obtain the cogeoidal (or geoidal, since the indirect effect is negligible) 
undulation N. 


Table 4. Computed Geo id undulations ( i2>o == 10*^ . m - 20 ) minus 

GEOS-3 Altimeter Undulations; Mean and RMS Differences 
in meters. 



Tonga Trench Area 

Indian Ocean Area 

Stdces' ftinction 

Mean 

Difference 

RMS 

Difference 

Mean 

Difference 

RMS 

Difference 

Unmodified 

0.8 

2.2 

0.8 

l.C 

Meissl's Mod. 

0.4 

1.1 

0.2 

0.7 

Molodenskii's mod. 
n = m = 20 

0.3 

0.9 

0.1 

0.6 

Molodenskli's mod. 
n = 10, m = 20 

0.4 

l.l 

0.3 

0.7 

i 


7. Colombo’s Method 


This section is devoted to the method developed by Colombo (1977), whose 
essential strategy is not to modify the error kernel in order to speed up the con- 
vergence rate of the error series, but to simply determine a kernel for the integra- 
tion over the cap which minimizes the truncation error. Under the condition of a 
"band-limited" gravity field (that is, one with a finite nmuber of harmonic compo- 
nents), the problem thus posed reduces to the solution of a finite number of unknown 
parameters from an equal number of linear equations. The (square) matrix to be 
inverted is extremely ill-conditioned, especially for the cap sizes of interest 
(^0 ^ 30°); nevertheless, Colombo (1977) achieved the desired result of significantly 
smaller truncation errors through a regularization of this matrix. In the following 
elaborations, analytic (recursive) expressions are derived for the decomposition of 
the matrix into lower and upper triangular matrices, leading further to analytic 
(recursive) expressions for the unknown quantities to be solved. The near singu- 
larity of the system is thereby not eliminated, but it will be shown that rigorously, 
without any form of regularization, the problem as stated above cannot be solved 
satisfactorily. In addition, if the gravity field is not "band-limited", then the error 
(without regularization) is far from minimal. 

The undulation is decomposed as follows 


N = 


R 

4ny 

R 

4TTy 


r SicoBiji) Ag da 

f 

Jkc(cos0) Agda+ j-j~^jAKc(cos^)Agda 
where by definition (the reason for this choice will be apparent below) 
Kc(y) = t ■^r^'^nPn(y) 


(126) 


(127) 


n-O 


SO that, evidently, we must have 


AKc (cosii)) = 


/ S(cos(£i) - Kc (cos(|)) , 
I S(cos0) , 


0 < 0 5 ^0 
00 < 0 ^ n 


(128) 


The unknown coefficients v„ are to be determined under the requirement of a 
minimum error. This error is given by the second integral in (126), which in view 
of (128) becomes 


6Nc 


" - j^JJkc(cos 0) Agda 

_R. f ^ f A 

" 2y ' 2y Ms 


(129) 

(130) 
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The coefficients 2nt„ are the eigenvalues of the first integral operator in (129) 
(see section 2), while the q„ are the Fourier coefficients of the kernel 


Rc(cos(£i) « |Kc(cos((i), 

0 0 ^ ^0 
^ n 

(131) 

that is, 



Qn ® r Kc(y) Pn(y) dy , 

0 s n < T5 

(132) 

Then clearly the eigenvalues of the second integral operator in (129) are 

2Trqr.. 


The result (130) is finally establ'shed by recalling our assumptions that Ago « 0 ® Agi . 
More compactly, we have 

R ^ 

“ 2 y ( tn “ Qn) ^gn (1^3) 

Considering (127) and (132), the coefficients q„ are functions of the parameters 
’ ir 

V’2r+1 __ . _ 

q« = — 2 — '''■ Jy 0 ^ n < n (134) 

For the moment, assume that the gravity field is composed of a finite number of 
harmonic components; i.e. there exists an integer n sue h that 

Ag„ = 0 , n > n (135) 

This bounding degree n is also used in the definition of Kc . Then from (133) 

R '''i 

fiNcs- = ~ li (tn-q„)Agn (136) 

^ A nans 

This is identically zero (the smallest possible error) if 

t„ = q„ , 0 5? n n (137) 

which with (134) becomes 

^ ^0 ~ ly Pn(y) dy , 0 s n s n (138) 

The (known) quantities t„ are given by (4). In matrix notation (138) is 

T = A V (139) 

where the positive definite matrix A (see Colombo (1977)) has elements 
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( 140 ) 


O n ^ 1 C ^ 

arn = R-(y) Pn(y) dy; r,n«=0 n 

*b 

and the vectors T, V contain elements t„ ( n » 0, . . , ,11 ) and Vr ( r » 0, . . . , n ) , 
respectively. Theoretically, under the condition (135), the system (139) of li -f 1 
equations can be solved for the TI <f 1 parameters Va , yielding a kernel Kc for 
which the truncation error 6Mcir ia zero. The gravity field in reality is described 
by an infinity of harmonics, but this only prevents the error from being absolutely 
zero. It is hoped, if (135) is ai^roximately true (say, if n » 180), that the 
actual erorr 

6Nc= -q„) Ag„ (141) 

remains small. 


Proceeding with the derivation of a recursive formula for the coefficients v^, 
we note that from (127) 


V 


n 



Kc(y) P„(y) dy 


(142) 


Kc can also be expanded in die interval [yo, 1] , where yo = cos (/iq ; this requires 
a corresponding orthogonal basis. Let 


y ~ £ z - Jl + I f where i = sin 


.3 


(143) 


then 


That is, the timctions P„((y + A - l)/;t) form an orthogonal basis for polynomials 
defined in [yo, 1]. Therefore, the expansion of Kc in [yo, 1) Is 

E 


Kc(y) = Kc(^z - £+ 1) 


n = 0 2 


Wn Pn( Z ) 


(145) 


This sum infinite since Kc is a polynomial of degree n in y , and hence in z . 
The Fourier coefficients are 

Wn = J ^ Kc ( - jg + 1) Pn( z ) dz , 0 ^ n s n (146) 

which upon inserting (127) become 


w„ = E 

r = 


2r+l 


0 2 
2r+l 


Vr Pp(^Z - A + 1) Pn( Z) dZ 


p = o 


grn Vr , 0 ^ n ^ n 


(147) 


where 
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(148) 


grn “ J ^ Pr(iz-A + l) Pb(«) dz ; r,n-0,...,n 

Since Pr( |z - jf, -f 1), when expanded in terms of P„(z), has no components beyond 
degree r , the orthogonality of Legendre polynomials guarantees that (cf. equation 


(108)) 



grn « 0 , n > r 


(149) 

Hence n 2r + l 

Wn “ L 9 Krn Vr , 

p=n * 

0 i n s: n 

(150) 

If we stipulate that 6NCH* 0 , 

then from (137) and (132) 


” I Kc(y) Pn(y) dy 
^0 




" Kc(Az - jt + 1) Pn(jtz - A + 1) dz (151) 

Substituting (145) yields 

I* ^ p 2 r + 1 

t„ = Jt o Wr Pr^Z) P„(Az-^+l)dz 

J-1 r*s 0 * 


P 2 r-f 1 

"" Hj Zij o gnr ^r» O^n^TT 

r a 0 ^ 


(152) 


where (149) has been used. Let W be a vector of dimension (n + 1) containing the 
coefficients w^ , and further let the (n + 1) x( n + 1) triangular matrices L and U 
be defined by 



In matrix notation, equations (152) and (150) can then be written concisely as 


T = LW, W = UV or T = LUV 


(155) 


With equation (139), we arrive at the desired result 
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A » LU 


(156) 


thus achieving a decomposition of matrix A into lower and upper triangular matrices. 


To determine the kernel Kc it is sufficient to solve only for W (see equation 
(145)). Inverting equation (152), it is easily verified that 


/2r + l . A . 2n-f 1 n 

Wr = ( 2 ^^'0 V’’ " ^ „5o 2 

Wo = to 


0 r n 


(157) 


Therefore, given values for ^ , r = 0, . . . ,TI (equation (4)) and a method to accurately 
compute g.„, equations (157) furnish (recursively) all coefficients w„. The following 
closed expression for g^ (cf. equation (113)) can be proved, e.g, , by mathematical 
induction: 


grn 


gfr 




+1 


0 < n <* r 


2jt‘~ 

2r + l ’ 


r S 0 


p = r-n-1 , q = r+n + 1 


(158) 

(159) 


However, because ( jt - 1)‘ is close to unity in magnitude (for the cap radii of 

Interest; e.g. (i>o = 10° -♦ .0076) and oscillates in sign, the summation in (158) 

on a digital computer is associated with a considerable loss in significant digits. 

A more suitable method to compute the g^„ in this case utilizes a recursion formula 
which is derived in Appendix B: 


gr+l,n 


Sr- l,n 


2r-fl 
2n + l 


g r, n + 1 ) » 


0 < n "S 


r + 1 




Sr ,0 = (Pr-l(yo) -Pr+l(yo)) , ^ > 0 

go,o = 2 ; gijO = 2(1-1) ; ~ 

grn = 0 , n > r 

Since ^<<1, grr (equation (159)) is an extremely small quantity for large r, while 
the values of gr,o are quite stable in magnitude as r increases; the values of gf^ 
vary between these bounds for 0 < n < r. Consequently, the coefficients w„ accor- 
ding to (157) are difficult to compute, being themselves large in magnitude and oscil- 
latory in sign. A considerable moderation in the extreme range of values of gp„ is 
effected by introducing a scaling factor Jt"" . Let 



Srn = gr„ i r,n ^ Q 


( 161 ) 


Then the recursive relationship (160) transforms into 
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K r + l,i» 


8r-i,n + 2n + l 8r, n + l) • 0 <" n S ffl ^ 


” ’ijlV+i) ‘ Pr+i(yo». r > 0 

Kcvo * 2 ; ^1,0 = 2(l-i) ; » j 

Kr,n " 0 , n> r 

Now (159) becomes "grr - 2/( 2 r+l) , r i: 0 , and with 
Wp » ;t"w„ , n 2: 0 
equation (157) changes simply to 

Wr = - t grn Wn , 0 < r <• n 

n“o " 

Wo = to 

The scaling described above does not alleviate totally the numerical difficulties 
since the values of g’m still vary over a broad range of magnitudes. For example, 
if (t>o = 10® , n = 36, then 

g88,o = -0. 119 { g36.i3 = 1. 82x10^^ ; g36.s4 = 5. 11x10^° ; gaa.aa = 0.0274 

and Wg = -1.17x10^'’} w^g = -4.20x10®; Wg 4 = 8.02x10^*5 Waa = 7. 51x10®^ 

The kernel Kc (equation (145)), in turn, acquires enormous oscillations, a charac- 
teristic that is both undesirable and detrimental to the actual truncation error. This 
error (equation (141)), as well as an indication of the computational accuracy of 
equations (162) and (164) are both provided through the evaluation of the coefficients 
q„. Substituting (143) and (145) into (142) results in 

qn = ^ r i Wr Pr( Z ) P„( ^ Z - i f 1) dz 

•^-1 r = 0 ^ 

= ^ Z Wr g„r , n S: 0 

r=sO ^ 

or, considering (161), (163), and (149) 

. f 2r+l 

q„ = jt L 5 gnr Wr , n & 0 

where M = min (n,lS) . For 0 ^ n s: If, we must have t„ =qn, where t„ is given 
by equation (4). Using the IBM System 370, model 168 with an AMDAHL 470 V16-II 
processor at The Ohio State University, Columbus, Ohio, with extended precision 
(more than 30 significant digits), the differences q„ - t„ were obtained with »{b = 10°, 
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(165) 

(166) 


>(162) 


(163) 


(164) 


W » 36 . Equal to zero for n « 0 , these differences deteriorated progressively until 
*aa “ 2. 3 X 10“® ; but then qg? - t*, 2. 2 x 10^°. The latter is essentially the 

first term in the truncation error. Therefore, while equations (137) are satisfied 
with at least seven digits of accuracy, the resulting strong oscillations of the kernel 
Kc are clearly reflected in the power spectrum of the error kernel, that is, the RMS 
truncation error, thus apparently ruling out any application of this method without 
some type of regularization of its ill-conditioned character. 


The determination of the coeffieients v„ 
since from (150) and (159) we have 


'n ^ fw„ - Y. 


2r+l 


g 


rn vp 


= - Y 


r s=n + l 
1 ! 


2r + l 


Ern Vr , 


r»n 4.1 


with 


Vr 




is numerically even more unpleasant, 

1 

0 ^ n is: n-l (167) 


Hence, these coefficients are excessively large (0o =* 10°, n = 5 |vn| is on the 
order of 10®® to 10®®), 


The rigorous solution for the coefficients v„ (or w„) as described above does 
not take into consideration the actual error 6 Nc (equation (141)) once the v„ are 
determined. Consider the RMS value of 6 Nc (equation (133)) s 



This quantity, which is to be minimized, resembles the square root of the sum of 
squared and weighted residuals In a typical least squares solution. To make the 
analogy strictly correct, the series in (168) is terminated at degree in > n . Then 
the coefficients to = ti » 0 , t„ = 2/(n - 1) , 2 s n m serve as the observations; 
the coefficients q„, 0 ^ n "s' m constitute the mathematical model ("adjusted 
observations"), being linear functions of the parameters v„, 0 ^ n ^ n. The 
degree variances c„ may be interpreted as weights. Note that the case m = n 
was discussed above. The degree of variability which enters the solution by allowing 
m > n may, through the minimization condition, have a mitigating influence on the 
magnitudes of the unknown parameters v„ . That is, the system of equations for the 
case m = n is known to be unstable, and the introduction of additional equations 
could have a regularizing effect on the system, since it is then overdetermined. 
However, the results of several computational experiments with m = 300, n = 50 
seem to indicate that the near singularity of the original system is not avertible 
with this approach. 
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8. Summary and ConcluBlon 


The geold undulation can be derived from a combination of gravity anomalies in 
a spherical cap (truncated Stokes' integral) and a (finite) set of harmonic coefficients. 
The resulting truncation error depends strongly on the kernel that is used for the 
integration. Several modifications to the kernel, as proposed by various authors, 
are scrutinized numerically in this report, These modifications are primarily 
designed to accelerate the convergence of the series which represents the error, 
thereby attempting to reduce its magnitude. It is found that these methods can claim 
such a reduction only if a number of low-degree harmonics (whose magnitudes are 
not directly controlled by the modifications of the kernel) can be deleted from the 
error and Incorporated in the computation of the undulation. This is not an absolute 
requirement since the error is also a function of the cap radius (lio (see Figure li); 
however, for small caps ( 5° < 0o 30®) , this seems to be the gene ml rule 
(Figures 6, 7, 9, 10). 

The various numerical investigations amply demonstrate that the integration 
of gravity anomalies in a cap of radius as small as 5® or 10®, supplemented by 
potential coefficients (e.g. GEM 9), can yield geoidal accuracies below the 1 meter 
level. For example, in Figure 10 (where integration and anomaly errors are neglected, 
but errors in GEM 9 are included), Meissl's (1971b) simple modification of the kernel 
gives an RMS error of 0.41 m for 10® i while a variation of Molodenskli's method 
(K= 10 ) produces an RMS error of 0. 33 m ( 0o ** 10°)> Geoid undulations, computed 
from a combination of actual gravity data in caps with t(io * 10® and the GEM 9 har- 
monic coefficients, have been compared to the corresponding GEOS-3 altimeter 
undulations which, in the areas considered, have an average accuracy of better than 
1 meter. While not directly verifying the results of Figure 10, this comparison 
does substantiate the significant (60%) improvement over the conventional method 
that can be achieved by suitable modifications of Stokes' function. It should be noted 
that the RMS errors of Figures 5, 6, 7, 9, and 10 are average global estimates that 
are based on the Tscheming and Rapp (1974) gravity model and a spherical approxi- 
mation of the geoid, the latter being no longer valid at the 20 cm to 30 cm level of 
accuracy. Furthermore, the errors of the gravity data, as well as errors associated 
with the numerical integration have been neglected entirely. 

Colombo's (1977) method of treating the truncation of Stokes' integral is included 
to show a different avenue of approach. Some numerical tests have been conducted, 
but this method requires further development and analysis. 
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If Off ■ J*"® Pr"(y) <*y • r a: 0 , than 

(2r+l)Crr - <2r-l)e,.i^.x+yolPr*(yo)+Rli<yo)J-2IV(yo)Pr-i(yo), r>0 (A.l) 

©0,0 ■ I + yo 

proof ! i 


The caee r 0 it eatily verified by noting that I{)(y ) • 1 • for all y . 
Suppoae that r > 0 , then with (30) we have 



rcrr » J^Pr(y) rp,(y) dy - J||°Pr(y)[y p;(y) -Pr’-i(y)l dy 

(A. 2) 

Let 

ar - J^°Pr(y) y Pr‘(y) dy . rio 

(A. 3) 

An integration by parts yields 



ar “ yPr®(y)|-l - J_||®Pr(y) tyPr'(y) +Pr(y)l dy 


Hence 

2 a, * yo P®( yo) + Pr®(-1) - e,. 

(A. 4) 

Putting (A. 3) into (A. 2), we obtain 



re„ - a, - J*^°Pr(y ) Pr’-i(y) dy 

(A. 5) 

Also 

(r+l)er+i,r+i » a,+x - J"° P,+x(y) Pr'(y)dy 



= a,^l - Pr+x(y)Pr(y)|i!l + f°Pr(y) PrVi(y)dy 

(A. 6) 


The last step follows from an integration by parts. Now add equations (A. 5) and 
(A. 6), and substitute (31): 

re„ + (r + l)e,^i^r+i = a, +ar+i - Pr+x<yo)Pr(yo) + Pr+i(-l)Pr(-l) + (2r+l) JV®(y) dy 
With (A. 4) and Pr(-l) = (-1)% this becomes 

(r+l+4)er+i^+i = (-r-i + 2r + l)err +iyoP®(yo) +^+^yoPr+i(yo) + 

- Pr+x(yo) Pr(Vo) - 1 
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which aimpUfiec to 

2 £ 12 o„,,,», - 2 !^e„ ♦lyo(ft’(yo)+p,!>(yo)l-Pr(yo)Pr«(yo) 
Equation (A. 1) followa immediately. 


Appendix B 

If - j Pr(i»-X+l) P„<a) d* , r,n s- 0 (B.X) 

then Rp+i,n * Kr-i,n + ^ ~ Kr»»+i)» {f > 0 (B- 2 ) 

"Pr+i<yo)l » r >0 (B.3) 

go,o “ 2 } gi,o *• 2(1-X) 5 gi^i “ f i (B.4) 

gr,n « 0 , n > r (B. 5) 

Proof; 


(B. 5) has already been eatablished in section 7, and equations (B.4) are easily 
verified. Using y ^ tz - i + 1 and (31) 

* 1 J “ A(2r-fi) ^^''-^<-^^^ " Pr+i(yo)l , r > 0 

% ' ' 

thus proving (B.3). 

r ^ 

Consider the integral j Pp(y) P»'(z) dz . integrating by parts, we obtain 
Pr^i(y)B.'(z)dz = 1 - Pr+i(-2«+l)(-l)"- iJ^PrXx(y)P„(z) dz (B. 6 ) 

Pr-i(y)Pn’(Z) dz = I -Pp_i(-2^ + l)(-l) - ^J_^Prli(y)Pn(Z) dz (B. 7) 

whexa the primes always denote differentiation with respect to the argument of the 
function. Subtracting (B.7) from (B. 6 ) and substituting (31) gives 

I ' tPr+l(y)-Pr-l(y)lPn’(Z)dZ = -Pr«(-2i + lK-l)"4-Pr-T(-2/ + lK-l)"-(2r + l) lgr„ 

(B. 8 ) 


- 49 - 


Also, U(31) is substituted lnto(B.l)i then 

!>’««(») - «'* 

Subtract {B. 10) from (B.9) and insert (B. 8), then 

-Kr-v,„ - 

- (-P(-2< + l)<-l)'>*' + lV.i(-2i+l)(-l)’’"^ - 
il2r + l), 

*“ 2n + l “ 8r»nnl 


(B.9) 
(B. 10) 

(2r+l)Ag,,„+i + 
(2r+l) igr,n-i)l 


thus proving also (B.2). 


